Ephedra regional gradient (ERG) analyses

ecoblender

alex.filazzola

Abstract

An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.

Hypothesis

We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.

Methods

load GAMs Gradient

Year<-c("2016","2016","2016","2016","2016","2016","2016","2017","2017","2017","2017","2017","2017","2017")
Site<-c("Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch","Barstow","Cuyama","HeartofMojave","PanocheHills","SheepholeValley","Tecopa","TejonRanch")
gamsq <-c(7992.36,6789.76,7903.21,5715.36,7779.24,7482.25,7447.69,6872.41,5882.89,7157.16,5595.04,7293.16,7106.49,6955.56)
gams.richard <-c(89.4,82.4,88.9,75.6,88.2,86.5,86.3,82.9,76.7,84.6,74.8,85.4,84.3,83.4)
gams <- c(89.6,71,84.7,62.7,88.3,80.6,80.7,63.6,51.4,73.4,49.5,69.4,68.3,65.8)

gams.data <- data.frame(Year,Site,gamsq,gams)

Weather for growing season

Climate patterns within study

season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))

## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)

## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)

### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days

season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)
##              Site frost.days
## 1         Barstow  24.725275
## 2          Cuyama  29.670330
## 3   HeartofMojave  25.824176
## 4    PanocheHills  10.439560
## 5 SheepholeValley   6.043956
## 6          Tecopa  23.626374
## 7      TejonRanch  50.549451
season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)
##              Site frost.days
## 1         Barstow  17.679558
## 2          Cuyama  19.889503
## 3   HeartofMojave  16.574586
## 4    PanocheHills  14.917127
## 5 SheepholeValley   1.785714
## 6          Tecopa  25.966851
## 7      TejonRanch  35.911602
## Both years had comparable number of frost days

## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}

season1.frost <- season1 %>% group_by(Site)  %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                  10
## 2          Cuyama                  10
## 3   HeartofMojave                  12
## 4    PanocheHills                   5
## 5 SheepholeValley                   7
## 6          Tecopa                  11
## 7      TejonRanch                  14
season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   7
## 4    PanocheHills                   6
## 5 SheepholeValley                   3
## 6          Tecopa                   7
## 7      TejonRanch                  13
## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  12.396694     5.750413
## 2          Cuyama  11.570248     3.352893
## 3   HeartofMojave  16.528926     4.542149
## 4    PanocheHills   2.479339     6.777686
## 5 SheepholeValley   2.479339     9.394167
## 6          Tecopa  11.570248     6.342149
## 7      TejonRanch  33.884298     1.438017
season2.frost <- season2 %>% group_by(Site) %>%  filter(year>2016) %>%  summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  11.666667     6.052500
## 2          Cuyama  16.666667     3.624167
## 3   HeartofMojave   8.333333     5.637838
## 4    PanocheHills   8.333333     6.065833
## 5 SheepholeValley   0.000000     9.386916
## 6          Tecopa  20.833333     5.938333
## 7      TejonRanch  28.333333     2.535833
season1.frost <- season1 %>% group_by(Site)  %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   4
## 4    PanocheHills                   2
## 5 SheepholeValley                   2
## 6          Tecopa                   4
## 7      TejonRanch                  10
season2.frost <- season2 %>% group_by(Site)  %>%  filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   3
## 4    PanocheHills                   3
## 5 SheepholeValley                   2
## 6          Tecopa                   5
## 7      TejonRanch                  10

Differences from Nutrient content

# nutrient data for each site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")


## join aridity with nutrients
nutrients.climate <- merge(nutrients,subset(gams.data, Year==2016), by=c("Site"))


## Nitrogen difference between shrub and sites
m.nit <- aov(log(N) ~ Site * microsite, data=nutrients)
summary(m.nit)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6  63.04   10.51  13.241 2.86e-09 ***
## microsite       1  40.16   40.16  50.614 2.24e-09 ***
## Site:microsite  6   4.96    0.83   1.042    0.408    
## Residuals      56  44.43    0.79                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.nit, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(N) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                diff      lwr      upr p adj
## shrub-open 1.514873 1.088317 1.941429     0
## Potassium difference between shrub and sites
m.pot <- aov(log(K) ~ Site * microsite, data=nutrients)
summary(m.pot)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 15.897   2.650  22.013 4.09e-13 ***
## microsite       1  4.445   4.445  36.934 1.14e-07 ***
## Site:microsite  6  1.605   0.268   2.223   0.0541 .  
## Residuals      56  6.740   0.120                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pot, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(K) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr       upr p adj
## shrub-open 0.5040081 0.3378755 0.6701406 1e-07
## Phosphorus difference between shrub and sites
m.pho <- aov(log(P) ~ Site * microsite, data=nutrients)
summary(m.pho)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 23.163   3.861   14.33 8.10e-10 ***
## microsite       1 10.546  10.546   39.15 5.79e-08 ***
## Site:microsite  6  3.394   0.566    2.10   0.0676 .  
## Residuals      56 15.085   0.269                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pho, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(P) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr      upr p adj
## shrub-open 0.7762734 0.5277317 1.024815 1e-07
nutrient.mean <- nutrients.climate %>% group_by(Site, gams) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)
 

ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=nitrogen), color="#E69F00", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#E69F00") + geom_jitter(aes(x=gams, y=potassium/20), color="#56B4E9", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=potassium/20), color="#56B4E9") + geom_jitter(aes(x=gams, y=phosphorus), color="#009E73", size=3) +ylab("soil nutrient content")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("rainfall continentality")+
annotate("segment", x = 72, xend = 74, y = 25, yend = 25,  colour = "#56B4E9", size=2) + ## add custom legend
  annotate("text", x = 65, y = 25,  label = "Nitrogen", size=5, hjust=0)+
annotate("segment", x = 72, xend = 74, y = 23, yend = 23,  colour = "#E69F00", size=2) + ## add custom legend
  annotate("text", x = 65, y = 23,  label = "Potassium", size=5, hjust=0)+
annotate("segment", x = 72, xend = 74, y = 21, yend = 21,  colour = "#009E73", size=2) + ## add custom legend
  annotate("text", x = 65, y = 21,  label = "Phosphorus", size=5, hjust=0)

m1 <- lm(nitrogen~poly(gams,2), data=nutrient.mean)
summary(m1)
## 
## Call:
## lm(formula = nitrogen ~ poly(gams, 2), data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  1.4326  3.1541 -1.6134 -1.2346  0.2699 -0.2336 -1.7750 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      5.4326     0.8325   6.526  0.00285 **
## poly(gams, 2)1  10.6460     2.2026   4.833  0.00844 **
## poly(gams, 2)2   9.0253     2.2026   4.098  0.01488 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.203 on 4 degrees of freedom
## Multiple R-squared:  0.9094, Adjusted R-squared:  0.8641 
## F-statistic: 20.08 on 2 and 4 DF,  p-value: 0.008208
m2 <- lm(phosphorus~gams, data=nutrient.mean)
summary(m2)
## 
## Call:
## lm(formula = phosphorus ~ gams, data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  3.9609  1.3929 -1.9844  0.5018 -0.6999 -5.4853  2.3140 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  23.6031    11.4769   2.057   0.0948 .
## gams         -0.1929     0.1432  -1.347   0.2357  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.399 on 5 degrees of freedom
## Multiple R-squared:  0.2664, Adjusted R-squared:  0.1196 
## F-statistic: 1.815 on 1 and 5 DF,  p-value: 0.2357
m3 <- lm(potassium~poly(gams,2), data=nutrient.mean)
summary(m3)
## 
## Call:
## lm(formula = potassium ~ poly(gams, 2), data = nutrient.mean)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##   8.0453 -13.5433   0.9712   4.7180 -15.7101  25.3691  -9.8502 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     186.557      6.705  27.824 9.93e-06 ***
## poly(gams, 2)1 -244.894     17.740 -13.805 0.000160 ***
## poly(gams, 2)2  159.200     17.740   8.974 0.000853 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.74 on 4 degrees of freedom
## Multiple R-squared:  0.9855, Adjusted R-squared:  0.9782 
## F-statistic: 135.6 on 2 and 4 DF,  p-value: 0.0002114

Native non-native comparisons

### annual plant community repsonse
spp.data  <- read.csv("Data/ERG.communitydata.csv")
spp.data[is.na(spp.data)] <- 0

## load community data
community <- read.csv("Data/ERG.communitydata.csv")
community[is.na(community)] <- 0


## load species list
spp.list <- read.csv("Data/ERG.specieslist.csv")

## convert data to long format
status <- gather(community, species, abundance, 13:53)
names(spp.list)[1] <- "species" ## rename species column for merge

## combine native-non-native status with data set
status <- merge(status, spp.list, by="species")

## summarize per plot native and non-natives, convert back to wide format
mean.status <- status %>% group_by(Year, Site, Microsite, Rep, status) %>% summarize(abundance=sum(abundance))
mean.status <- spread(mean.status, status, abundance)
status.data <- mean.status

## calculate site means for native and non-native
mean.status <- mean.status %>% group_by(Year, Site, Microsite) %>% summarize(native=mean(native), non.native=mean(non.native))

Phylogenetics

library(picante)
library(ape)
library(brranching)
library(taxize)

## load species list to create tree
spp.list <- read.csv("Data/ERG.specieslist.csv")

## create a combined species name column and drop erinoginum because not to species level
taxa <- paste(spp.list$Genus, spp.list$Species.name)
taxa <- taxa[taxa!="Erinoginum spp"]

## create tree of phylogeny
tree <- phylomatic(taxa=taxa, get = 'POST')

## calculate branch lengths using  Grafen's method
## Grafen, A. (1989) The phylogenetic regression. Philosophical Transactions of the Royal society of London. Series B. Biological Sciences, 326, 119–157.
tree2 <- compute.brlen(tree, method="Grafen")

## view tree and outputs to verify branches
plot(tree2)

tree2
## 
## Phylogenetic tree with 39 tips and 29 internal nodes.
## 
## Tip labels:
##  cryptantha_barbigera, cryptantha_intermedia, phacelia_crenulata, phacelia_tanacetifolia, pectocarya_spp, amsinckia_grandiflora, ...
## Node labels:
##  poales_to_asterales, eudicots, , , asterids, ericales_to_asterales, ...
## 
## Rooted; includes branch lengths.
## format community to only have site and microsite
comm <- community %>% group_by(Year, Site, Microsite) %>% summarise_if(is.numeric, funs(sum(., na.rm=T)))
comm <- comm[,c(1:3,13:53)] ## extract species abundances only
comm <- data.frame(comm) ## convert to dataframe
comm[,"site.micro"] <- paste(comm$Site,comm$Microsite,as.character(comm$Year)) ## create site by microsite column
rownames(comm) <- comm[,length(comm)] ## replace row names with site by microsite
comm[,c("Year","Site","Microsite","Erinoginum.spp","Boraginaceae.sp.","site.micro")] <- NULL ## drop all columns but ones for analysis

## match species name format
spp.list[,"spp.name"] <- paste(spp.list$Genus, spp.list$Species.name)
colnames(comm) <- spp.list$spp.name[match(colnames(comm), spp.list$Species.shorthand)]

## match species format
names(comm) <- gsub(" ", "_", names(comm), fixed=T) ## underscores instead of spaces
names(comm) <- gsub("__", "_", names(comm), fixed=T) ## replace double underscores with one
names(comm) <- tolower(names(comm)) ## lower case names


## mean phylogenetic distance
mpd.data <- mpd(comm, cophenetic(tree2), abundance.weighted=T)

## format dataframe
mpd.data <- data.frame(Year=c(rep(2016,14),rep(2017,14)),Microsite=c("open","shrub"),site=rownames(comm), mpd=mpd.data)
mpd.data[,"Site"] <- gsub(" open", "", mpd.data[,"site"]) 
mpd.data[,"Site"] <- gsub(" shrub", "", mpd.data[,"Site"]) 
mpd.data[,"Site"] <- gsub(" 2016", "", mpd.data[,"Site"]) 
mpd.data[,"Site"] <- gsub(" 2017", "", mpd.data[,"Site"]) 
mpd.data[,"site"] <- NULL
mpd.data[is.na(mpd.data)] <- 0 ## barstow no plants
 

## getspecies richness for sites
spp.rich <- community %>% group_by(Year, Site, Microsite) %>% summarize(richness=mean(Richness))
spp.rich <- data.frame(spp.rich)

## join data
mpd.rich <- merge(mpd.data, spp.rich, by=c("Site","Microsite","Year"))

Community analysis

### indicator species analysis
library(indicspecies)

## indicator species analysis
indval <- multipatt(community[,13:53], community$Microsite, control=how(nperm=999))
summary(indval)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: IndVal.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 41
##  Selected number of species: 14 
##  Number of species associated to 1 group: 14 
## 
##  List of species associated to each combination: 
## 
##  Group open  #sps.  8 
##                   stat p.value    
## E.cicutarium     0.558   0.001 ***
## L.gracilis       0.427   0.001 ***
## A.wrangelianus   0.201   0.001 ***
## C.brevipes       0.171   0.032 *  
## A.lentiginosus   0.149   0.027 *  
## C.rigida         0.149   0.023 *  
## A.grandiflora    0.129   0.018 *  
## Boraginaceae.sp. 0.120   0.020 *  
## 
##  Group shrub  #sps.  6 
##                  stat p.value    
## C.fremontii     0.379   0.001 ***
## M.glabrata      0.277   0.001 ***
## P.tanacetifolia 0.267   0.001 ***
## A.tessellata    0.199   0.005 ** 
## B.nigra         0.179   0.014 *  
## E.nitidum       0.172   0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
site.names <- as.character(unique(community$Site))
community.site <- subset(community, Site==site.names[7])



## species accumilation curves
accum1 <- specaccum(subset(community.site, Microsite=="shrub", select=13:53), method="random",  permutations=1000)
accum2 <- specaccum(subset(community.site, Microsite=="open", select=13:53), method="random",  permutations=1000)


par(mar=c(4.5,4.5,0.5,0.5))
plot(accum2,  col="#E69F00",ci.col="#E69F0090", xlab="Sampling effort", cex.axis=1.5, cex.lab=1.8, ylab="Species richness", lwd=2, ylim=c(0,16))
plot(accum1,add=T , col="#56B4E9",   ci.col="#56B4E990", lwd=2)

micros <- community$Microsite[1:120]

indval <- multipatt(community.site[,13:53], micros, control=how(nperm=999))
summary(indval)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: IndVal.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 41
##  Selected number of species: 4 
##  Number of species associated to 1 group: 4 
## 
##  List of species associated to each combination: 
## 
##  Group open  #sps.  1 
##                   stat p.value  
## Boraginaceae.sp. 0.316   0.029 *
## 
##  Group shrub  #sps.  3 
##               stat p.value    
## A.tessellata 0.531   0.001 ***
## B.nigra      0.475   0.013 *  
## M.affinis    0.343   0.048 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cluster analysis

library("cluster")
library("factoextra")
library("magrittr")

comm.sum <- community %>%  group_by(Microsite, Site) %>% summarize_if(is.numeric, funs(sum))
comm.sum <- data.frame(comm.sum)
rownames(comm.sum) <- paste(comm.sum$Microsite,comm.sum$Site)

## coreelation matrix

### Distance measures
res.dist <- get_dist(comm.sum[,13:53], stand = TRUE, method = "euclidean")
fviz_dist(res.dist,    gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

## dendrogram

## transform data
comm.trans <- decostand(comm.sum[,13:53], "hell")

## clean data for ordination
## see distribution of spp
boxplot(comm.trans, xaxt="n")
labs <- colnames(comm.trans)
text(cex=0.8, x=1:41-1, y=-0.12, labs, xpd=TRUE, srt=45)

## remove spp with only one instancve
comm.trans <- comm.trans[,!colSums(comm.trans)==apply(comm.trans, 2, max)]


## replace outliers with mean
avg.max <- function(x) {
  y =  max(x)
  avg = mean(x)
  x[x==y] <- avg
  return(x)
}
  
## dot chart to find outliers (remove or Avg)
#dotchart(comm.trans$E.cicutarium) ## Erodium
comm.trans[,"E.cicutarium"] <- avg.max(comm.trans$E.cicutarium) ## removed outlier
#dotchart(comm.trans$A.wrangelianus) ## Chilean trefoil
comm.trans[,"A.wrangelianus"] <- NULL ## removed entire species
#dotchart(comm.trans$A.lentiginosus) ##
comm.trans[,"A.lentiginosus"] <- avg.max(comm.trans$A.lentiginosus) ## removed outlier
#dotchart(comm.trans$H.vulgare) ## Hordeum
comm.trans[,"H.vulgare"] <- NULL ## removed entire species
#dotchart(comm.trans$Pectocarya.spp) ## 
comm.trans[,"Pectocarya.spp"] <- avg.max(comm.trans$Pectocarya.spp) ## removed outlier
#dotchart(comm.trans$L.arizonicus) ## Lupin 
comm.trans[,"L.arizonicus"] <- NULL ## removed entire species
#dotchart(comm.trans$M.bellioides) ##
comm.trans[,"M.bellioides"] <- NULL ## removed entire species
#dotchart(comm.trans$B.nigra) ## Mustard
comm.trans[,"B.nigra"] <- NULL ## removed entire species
#dotchart(comm.trans$N.demissum) ## Sand Matt
comm.trans[,"N.demissum"] <- avg.max(comm.trans$N.demissum) ## removed outlier
#dotchart(comm.trans$L.gracilis) ## 
comm.trans[,"L.gracilis"] <- avg.max(comm.trans$L.gracilis) ## removed outlier
#dotchart(comm.trans$A.grandiflora) ##  Amsinckia
comm.trans[,"A.grandiflora"] <- NULL ## removed entire species
#dotchart(comm.trans$B.diandrus) ##  Bromus diandrus
comm.trans[,"B.diandrus"] <- NULL ## removed entire species
#dotchart(comm.trans$M.affinis) ##  
comm.trans[,"M.affinis"] <- NULL ## removed entire species
#dotchart(comm.trans$P.crenulata) ##  
comm.trans[,"P.crenulata"] <- NULL ## removed entire species
#dotchart(comm.trans$Erinoginum.spp) # buckwheat
comm.trans[,"Erinoginum.spp"] <- avg.max(comm.trans$Erinoginum.spp) ## removed outlier
#dotchart(comm.trans$C.lasiophyllus) # 
comm.trans[,"C.lasiophyllus"] <- NULL ## removed entire species



library(usdm)
library(corrplot)

## Check for collinearity
vifcor(comm.trans)
## 4 variables from the 25 input variables have collinearity problem: 
##  
## E.wallacei P.insularis E.glyptosperma L.decorus 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( L.gracilis ~ S.barbatus ):  -0.003266766 
## max correlation ( C.fremontii ~ M.glabrata ):  0.8450794 
## 
## ---------- VIFs of the remained variables -------- 
##          Variables VIF
## 1         B.rubens Inf
## 2     E.cicutarium Inf
## 3       S.barbatus Inf
## 4      C.barbigera Inf
## 5    C.claviformis Inf
## 6   A.lentiginosus Inf
## 7         C.rigida Inf
## 8  P.tanacetifolia Inf
## 9   Pectocarya.spp Inf
## 10     D.capitatum Inf
## 11      C.brevipes Inf
## 12      M.glabrata Inf
## 13     C.fremontii Inf
## 14    C.intermedia Inf
## 15      N.demissum Inf
## 16      L.gracilis Inf
## 17    A.tessellata Inf
## 18       P.secunda Inf
## 19  L.sparsiflorus Inf
## 20  Erinoginum.spp Inf
## 21       E.nitidum Inf
## drop Claviformis because of collinear problems in P.insularis
comm.trans[,"C.claviformis"] <- NULL ## removed entire species
comm.trans[,"L.decorus"] <- NULL ## removed entire species
comm.trans[,"E.glyptosperma"] <- NULL ## removed entire species
comm.trans[,"E.wallacei"] <- NULL ## removed entire species

cor.dat <- cor(comm.trans[,1:20])
corrplot(cor.dat, method="number")

## calculate distances
dis <- vegdist(comm.trans, method="bray")


## cluster analysis
clus <- hclust(dis, "ward.D2")

## cluster colours
colz <- c("#3e4444", "#86af49", "#405d27","#c1946a")

fviz_dend(clus, k = 4, # Cut in four groups
          cex = 0.7, cex.axis=2, # label size
          k_colors = colz,
          color_labels_by_k = TRUE, # color labels by groups
          rect = TRUE # Add rectangle around groups
          ) 

grp <- cutree(clus, 4)

## Interpretation of classes

## get soil moisture from site/microsite

## Soil moisture
smc <- read.csv("Data/ERG.phytometer.census.csv")
smc.end <- subset(smc, Census=="end")
smc.avg <- smc.end %>% group_by(Site, Microsite) %>% summarize(mean.smc=mean(swc, na.rm=T))
smc.avg <- data.frame(smc.avg)
smc.avg[,"grp"] <- as.factor(c(3,3,2,1,4,4,1,1,4,4,3,3,2,2))

boxplot(mean.smc ~ grp, data=smc.avg)

## CA or PCA
dca1 <- decorana(comm.trans) ## length of gradient >2 & determine relative differences in community composition 


## conduct ordination
ord <- cca(comm.trans)
summary(ord)
## 
## Call:
## cca(X = comm.trans) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total            1.53          1
## Unconstrained    1.53          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.6615 0.2873 0.1656 0.10566 0.08765 0.07595 0.05451
## Proportion Explained  0.4324 0.1878 0.1082 0.06907 0.05729 0.04964 0.03563
## Cumulative Proportion 0.4324 0.6202 0.7284 0.79749 0.85478 0.90443 0.94006
##                           CA8     CA9    CA10     CA11     CA12     CA13
## Eigenvalue            0.03070 0.02546 0.01821 0.008974 0.004570 0.003798
## Proportion Explained  0.02006 0.01664 0.01190 0.005866 0.002987 0.002482
## Cumulative Proportion 0.96012 0.97676 0.98866 0.994530 0.997518 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                     CA1      CA2      CA3      CA4      CA5       CA6
## B.rubens        -1.5272  0.37620  0.03693 -0.17971  0.26195 -0.007680
## E.cicutarium    -0.5874 -0.19213 -0.20338 -0.20206 -0.36108  0.524976
## S.barbatus       0.1974 -0.64533  0.21879  0.05816  0.09674 -0.053334
## C.barbigera      0.4900 -0.11013  0.11027  0.09640  0.08716  0.152044
## A.lentiginosus   0.2199 -0.28948 -0.59265 -0.08500  0.77156 -0.296433
## C.rigida         0.4907 -0.59656 -0.53938  0.07713  0.37931 -0.042119
## P.tanacetifolia -1.2902  0.91318 -0.45764  2.53829  0.03756 -0.152588
## Pectocarya.spp   0.5495 -0.11384 -0.59153 -0.11011  0.04276  0.166931
## D.capitatum     -0.7173 -0.38053 -0.34301  0.92059 -1.47031  1.938948
## C.brevipes       0.4996 -0.08694 -1.06985 -0.06884 -0.31678  0.438663
## P.insularis      0.8031  0.56576 -0.62483 -0.12389  0.05743 -0.096057
## M.glabrata       0.7884  0.84798  0.94548 -0.09754 -0.18919  0.009159
## C.fremontii      0.8264  0.71519  0.19154 -0.07910 -0.11619  0.016657
## C.intermedia     0.3301 -0.17327 -0.89534 -0.13610 -0.03923 -0.580017
## N.demissum       0.7716  0.44620 -0.72001  0.09033  0.21708 -0.241934
## L.gracilis      -0.9869 -0.31033  0.07421  0.35879 -1.02622 -0.476879
## A.tessellata    -1.0035 -0.36194  0.26967 -0.84774 -1.13864 -1.238336
## P.secunda       -1.5685  0.52156  0.01749  0.80285  0.12080 -0.361127
## L.sparsiflorus   0.7841  1.12673  0.93614  0.24022  0.04660  0.069996
## Erinoginum.spp   0.5728 -1.00733  0.32438  0.36758  1.12435 -0.440374
## E.nitidum        0.8450  0.85547 -0.32831 -0.02364 -0.06814 -0.220852
## 
## 
## Site scores (weighted averages of species scores)
## 
##                           CA1     CA2      CA3      CA4      CA5       CA6
## open Barstow           0.3623 -1.3253  0.59499  0.23046  0.56897  0.613534
## open Cuyama           -0.2820 -0.9387 -0.39810  0.01254 -2.11537  2.774795
## open HeartofMojave     0.7431  0.2335 -2.05624 -0.52149  0.05849  0.033486
## open PanocheHills     -1.6700  0.4724  0.08096 -1.23638  1.43764  0.474226
## open SheepholeValley   0.9186  0.9037  0.77143  0.17618  0.30513  0.116582
## open Tecopa            0.5210 -1.3909 -0.13250  0.41626  1.33826 -0.340390
## open TejonRanch       -0.6455 -0.9594  0.23131 -0.09915 -1.38186 -1.124251
## shrub Barstow          0.5555 -0.5840  0.87451  0.02079  0.01283  0.724468
## shrub Cuyama          -1.7597  0.9562 -0.21109  3.09525  0.07451 -0.062796
## shrub HeartofMojave    0.7511  0.7261 -1.52964  0.11670 -0.12309 -0.542992
## shrub PanocheHills    -2.0371  0.9309 -0.05454 -1.74126  1.62903  1.240252
## shrub SheepholeValley  1.0489  1.6807  1.57641 -0.25781 -0.55721 -0.007623
## shrub Tecopa           0.4692 -1.4952  0.52997  0.40750  1.29037 -0.766954
## shrub TejonRanch      -1.1319 -0.1862  0.35145 -1.12008 -1.13771 -1.730813
## calculate priority
spp.priority <- colSums(comm.trans)

## plot ordination
par(mar=c(4.5,4.5,0.5,0.5))
plot(ord, type="n")
ordiellipse(ord, grp, lty = 2, col = "grey80", draw="polygon", alpha=150)
orditorp(ord, display = "species", cex = 0.7, col = "darkorange3", priority=spp.priority, air=0.5)
orditorp(ord, display = "sites", cex = 0.7, col = "darkslateblue", air=0.1)


##collect environmental variables for site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
nutrients.mean <- nutrients %>% group_by(microsite, Site) %>% summarise_all(funs(mean))
nutrients.vars <- data.frame(nutrients.mean)


## summarize daily means across sites
HOBO <- read.csv("Data/ERG.logger.data.csv")
HOBO.data <- HOBO %>%  group_by(Site, Microsite, Year, Month, Day) %>% summarize(Temp=mean(Temp),RH=mean(RH))

means <- HOBO.data %>% group_by(Microsite, Site) %>% summarize(Temp.=mean(Temp, na.rm=T),rh=mean(RH, na.rm=T),temp.se=se(Temp),rh.se=se(RH))
means <- data.frame(means)

## soil moisture
smc.early <- subset(smc, Census=="emergence")
smc.avg <- smc.early %>% group_by(Microsite, Site) %>% summarize(SMC=mean(swc, na.rm=T))

envs <- data.frame(swc=smc.avg[,"SMC"],nutrients.vars[,c("N","P","K")], means[,c("Temp.","rh")])

## Check for collinearity
cor(envs)
##              SMC           N           P           K      Temp.         rh
## SMC    1.0000000 -0.31692385  0.44854712  0.81694831 -0.6372354  0.9185204
## N     -0.3169238  1.00000000  0.05879919 -0.05275422  0.4167329 -0.4367433
## P      0.4485471  0.05879919  1.00000000  0.64885814 -0.5303896  0.5342666
## K      0.8169483 -0.05275422  0.64885814  1.00000000 -0.5769100  0.7750256
## Temp. -0.6372354  0.41673289 -0.53038959 -0.57690998  1.0000000 -0.8650006
## rh     0.9185204 -0.43674328  0.53426662  0.77502564 -0.8650006  1.0000000
envs[,"K"] <- NULL ## drop potassium for being correlated
envs[,"rh"] <- NULL ## drop humidity for being correlated

ord.env <- envfit(ord, envs)
plot(ord.env)

Ordination including aridity and microsite (pCCA)

Partial Constrained Correspondence Analysis for community data

Compare Ambient Community - mixed model results

community.arid <- merge(community, gams.data, by=c("Year","Site"))

mean.spp <- community.arid %>% group_by(Microsite, Site, gams, Year) %>%  summarize(abd=mean(Abundance),rich=mean(Richness), bio=mean(Biomass))

mean.status2 <- merge(mean.status,gams.data, by=c("Year","Site"))

m1 <- lmer(ihs(bio)~Microsite * gams + (1|Year), data=mean.spp)
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.9789, p-value = 0.8234
anova(m1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Microsite       1.6259  1.6259     1 23.000  11.7820  0.002272 ** 
## gams           18.1519 18.1519     1 23.367 131.5327 4.448e-11 ***
## Microsite:gams  1.0755  1.0755     1 23.000   7.7937  0.010365 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##      R2m      R2c 
## 0.549095 0.945669
m2 <- lmer(ihs(abd)~Microsite * gams + (1|Year), data=mean.spp)
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.97267, p-value = 0.6537
car::Anova(m2, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ihs(abd)
##                  Chisq Df Pr(>Chisq)    
## Microsite       0.1875  1     0.6650    
## gams           63.2057  1  1.862e-15 ***
## Microsite:gams  0.1083  1     0.7421    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m2.nat <- lmer(ihs(native)~Microsite * poly(gams,2) + (1|Year), data=mean.status2)
shapiro.test(resid(m2.nat))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2.nat)
## W = 0.98084, p-value = 0.8701
anova(m2.nat)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Microsite               1.2604  1.2604     1    22  1.1237 0.30063  
## poly(gams, 2)           8.0514  4.0257     2    22  3.5890 0.04477 *
## Microsite:poly(gams, 2) 2.5294  1.2647     2    22  1.1275 0.34184  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m2.exo <- lmer(ihs(non.native)~Microsite * gams+ (1|Year), data=mean.status2)
shapiro.test(resid(m2.exo))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2.exo)
## W = 0.90696, p-value = 0.01672
anova(m2.exo)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Microsite       0.247   0.247     1 23.000  0.1831    0.6727    
## gams           74.389  74.389     1 23.965 55.0565 1.175e-07 ***
## Microsite:gams  0.175   0.175     1 23.000  0.1294    0.7223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.6343625 0.8521801
m3 <- lmer(rich~Microsite * poly(gams,2) + (1|Year), data=mean.spp)
shapiro.test(resid(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m3)
## W = 0.93835, p-value = 0.1004
car::Anova(m3, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: rich
##                           Chisq Df Pr(>Chisq)   
## Microsite                0.6282  1   0.428012   
## poly(gams, 2)           10.9223  2   0.004249 **
## Microsite:poly(gams, 2)  2.5197  2   0.283699   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values 
##       R2m       R2c 
## 0.3445392 0.3495475
## Merge MPD with gams
mpd.gams <- merge(mpd.rich, gams.data, by=c("Site","Year"))

m4 <- lmer(mpd~Microsite * gams + (1|Year), data=mpd.gams)
shapiro.test(resid(m4))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m4)
## W = 0.98022, p-value = 0.8555
anova(m4)
## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Microsite      0.67216 0.67216     1    24  5.5363 0.02716 *
## gams           0.00313 0.00313     1    24  0.0258 0.87369  
## Microsite:gams 0.54576 0.54576     1    24  4.4952 0.04452 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m4) ## R squared values 
##       R2m       R2c 
## 0.2060199 0.2060199
### Biomass
p1 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=ihs(bio)), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(bio)), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#E69F00", fill="#E69F0080", data=subset(mean.spp, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,3),aes(x=gams, y=ihs(bio)), color="#56B4E9", fill="#56B4E970", data=subset(mean.spp, Microsite=="shrub"), lwd=1.4)  + ylab("Annnual plant biomass")+ xlab("")  #annotate("text", x=-4.5,y=4, label="a", size=8)+coord_cartesian(ylim=c(0, 4))

### abundance
p2 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=ihs(abd)), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(abd)), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open"))  + theme_Publication()+  stat_smooth(method="lm", formula= ihs(y)~x,aes(x=gams, y=abd), color="black", fill="Grey50", data=mean.spp, lwd=1.4)+ ylab("Annual plant abundance")+ xlab("")+ annotate("text", x=-4.5,y=6, label="b", size=8)

### richness
p3 <- ggplot(mean.spp) + geom_jitter(aes(x=gams, y=rich), color="#56B4E9", size=3, data=subset(mean.spp, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=rich), color="#E69F00", size=3, data=subset(mean.spp, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#E69F00", fill="#E69F0080", data=subset(mean.spp, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=rich), color="#56B4E9", fill="#56B4E970", data=subset(mean.spp, Microsite=="shrub"), lwd=1.4)  + ylab("Species richness") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=4, label="c", size=8)

### phylogenetic diveristy
p4 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#56B4E9", size=3, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#E69F00", size=3, data=subset(mpd.gams, Microsite=="open"))  +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#56B4E9",fill="#56B4E980", lwd=1.4)+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#E69F00",fill="#E69F0080", lwd=1.4)+ ylab("Phylogenetic Community Dissimilarity") + xlab("Aridity Gradient")+ theme_Publication()+ annotate("text", x=-4.5,y=1.8, label="d", size=8)

grid.arrange(p1,p2,p3,p4)

Compare environmental differences - mixed model results

### Nutrients 
nutrient.mean <- nutrients.climate %>% group_by(Site, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))

## potassium
m5 <- lm(potassium ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)
## 
## Call:
## lm(formula = potassium ~ poly(gams, 2) * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.224 -15.052  -5.456   7.166  51.901 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     141.46      12.13  11.662 2.67e-06 ***
## poly(gams, 2)1                 -217.98      45.39  -4.803 0.001351 ** 
## poly(gams, 2)2                  150.46      45.39   3.315 0.010617 *  
## micrositeshrub                   90.20      17.15   5.258 0.000766 ***
## poly(gams, 2)1:micrositeshrub  -256.71      64.19  -4.000 0.003952 ** 
## poly(gams, 2)2:micrositeshrub   149.36      64.19   2.327 0.048389 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.09 on 8 degrees of freedom
## Multiple R-squared:  0.9641, Adjusted R-squared:  0.9416 
## F-statistic: 42.95 on 5 and 8 DF,  p-value: 1.438e-05
anova(m5)
## Analysis of Variance Table
## 
## Response: potassium
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## poly(gams, 2)            2 170636   85318  82.835 4.502e-06 ***
## microsite                1  28476   28476  27.648 0.0007663 ***
## poly(gams, 2):microsite  2  22053   11026  10.706 0.0054741 ** 
## Residuals                8   8240    1030                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m5) ## R squared values 
##       R2m       R2c 
## 0.9429147 0.9429147
p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=potassium), color="#56B4E9", size=3, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=potassium), color="#E69F00", size=3, data=subset(nutrient.mean, microsite=="open"))+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#56B4E9", fill="#56B4E980", data=subset(nutrient.mean, microsite=="shrub"), lwd=1.4) +  
  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=potassium), color="#E69F00", fill="#E69F0080",data=subset(nutrient.mean, microsite=="open"), lwd=1.4)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ annotate("text", x=-4.5,y=32, label="a", size=8)+coord_cartesian(ylim=c(0, 32))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
se <- function(x) sd(na.omit(x))/sqrt(length(na.omit(x)))

swc <- smc %>% group_by(Year, Site, Microsite, Census) %>% summarize(swc.avg=mean(swc), swc.se=se(swc))
smc.gams <- merge(swc, gams.data, by=c("Site","Year"))


## early swc
smc.early <- subset(smc.gams, Census=="emergence")


m6 <- lmer(swc.avg~Microsite * poly(gams,2) + (1|Year), data=smc.early)
shapiro.test(resid(m6))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m6)
## W = 0.91432, p-value = 0.02516
anova(m6)
## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Microsite                  1.53    1.53     1 21.000  0.0837    0.7752    
## poly(gams, 2)           1449.15  724.58     2 21.343 39.5604 6.617e-08 ***
## Microsite:poly(gams, 2)    3.69    1.85     2 21.000  0.1007    0.9046    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m6) ## R squared values 
##       R2m       R2c 
## 0.5464095 0.8953768
### soil moisture at emergence
p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#56B4E9", size=3, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#E69F00", size=3, data=subset(smc.early, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=1.4, color="black") + ylab("Soil moisture at emergence (%)") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=30, label="c", size=8)
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## Create seasons for hobo
season1 <- subset(HOBO.data, Year==2015 & Month > 10 | Year==2016 & Month < 5)
season2 <- subset(HOBO.data, Year==2016 & Month > 10 | Year==2017 & Month < 5)
HOBO.season <- rbind(season1,season2)
HOBO.season[,"season"] <- c(rep("season.1",nrow(season1)),rep("season.2",nrow(season2)))
HOBO.season <- na.omit(HOBO.season)

hobo.means <- HOBO.season %>% group_by(season, Microsite, Site) %>% summarize(temp.avg=mean(Temp),temp.var=var(Temp),temp.sd=sd(Temp),temp.se=se(Temp),rh.avg=mean(RH, na.rm=T),rh.var=var(RH, na.rm=T),rh.se=se(RH),frost=(sum(na.omit(Temp)<0)/length(na.omit(Temp))*100),heat=(sum(na.omit(Temp)>30)/length(na.omit(Temp))*100))
hobo.means[,"Year"] <- c(rep(2016,14),rep(2017,14))

hobo.arid <- merge(hobo.means, gams.data, by=c("Site","Year"))
# hobo.arid[,"CV"] <-hobo.arid[,"temp.sd"] /hobo.arid[,"CV"] 

## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m7)
## W = 0.97826, p-value = 0.8069
anova(m7)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Microsite                 1.51    1.51     1 21.000  0.2547    0.6190    
## poly(gams, 2)           707.64  353.82     2 21.104 59.7042 2.051e-09 ***
## Microsite:poly(gams, 2)   0.11    0.05     2 21.000  0.0092    0.9908    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m7) ## R squared values 
##       R2m       R2c 
## 0.3627874 0.9560767
p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#56B4E9", size=3, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#E69F00", size=3, data=subset(hobo.arid, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#E69F00", fill="#E69F0080", data=subset(hobo.arid, Microsite=="open"), lwd=1.4) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#56B4E9", fill="#56B4E970", data=subset(hobo.arid, Microsite=="shrub"), lwd=1.4)  + ylab("Temperature variation") + xlab("Aridity Gradient")+ annotate("text", x=-4.5,y=125, label="d", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
shrubz <- read.csv("Data/ERG.shrub.csv")
shrub.mean <- shrubz %>% group_by(Site, Microsite) %>% summarize(compact=mean(Compaction))

shrub.clim <- data.frame(shrub.mean, gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))
#shrub.clim <- subset(shrub.clim, Year==2016)

m8 <- lm(log(compact)  ~  Microsite * gams, data=shrub.clim)
summary(m8)
## 
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37348 -0.12596  0.04006  0.12843  0.33669 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.138777   0.804662  -0.172   0.8665  
## Micrositeshrub      -3.101869   1.137963  -2.726   0.0213 *
## gams                 0.007791   0.010038   0.776   0.4556  
## Micrositeshrub:gams  0.036028   0.014196   2.538   0.0295 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared:  0.6967, Adjusted R-squared:  0.6058 
## F-statistic: 7.658 on 3 and 10 DF,  p-value: 0.006
anova(m8)
## Analysis of Variance Table
## 
## Response: log(compact)
##                Df  Sum Sq Mean Sq F value   Pr(>F)   
## Microsite       1 0.18836 0.18836  3.3176 0.098548 . 
## gams            1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams  1 0.36569 0.36569  6.4409 0.029469 * 
## Residuals      10 0.56776 0.05678                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m8) ## R squared values 
##       R2m       R2c 
## 0.6386405 0.6386405
p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=compact), color="#56B4E9", size=3, data=subset(shrub.clim, Microsite=="shrub"))+ geom_jitter(aes(x=gams, y=compact), color="#E69F00", size=3, data=subset(shrub.clim, Microsite=="open"))+  stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#56B4E9", fill="#56B4E980", data=subset(shrub.clim, Microsite=="shrub"), lwd=1.4) +  
  stat_smooth(method="lm", formula= y~exp(x),aes(x=gams, y=compact), color="#E69F00", fill="#E69F0080",data=subset(shrub.clim, Microsite=="open"), lwd=1.4)+ ylab("Soil compaction") + xlab("")+ theme_Publication() + annotate("text", x=-4.4,y=2.5, label="b", size=8) +coord_cartesian(ylim=c(0.5, 2.5))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p5,p8,p6,p7)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

### Microsite comparisons



# nut <- read.csv("Data//ERG.soilnutrients.csv")
# 
# fit1 <- aovp(N ~ microsite, data=nut)
# summary(fit1)
# t.test(N ~ microsite, data=nut)
# fit2 <- aovp(P ~ microsite, data=nut)
# t.test(P ~ microsite, data=nut)
# summary(fit2)
# fit3 <- aovp(K ~ microsite, data=nut)
# t.test(K ~ microsite, data=nut)
# summary(fit3)
# 
# compact <- read.csv("Data//ERG.shrub.csv")
# 
# fit4 <- aovp(Compaction ~ Microsite, data=compact)
# t.test(Compaction ~ Microsite, data=compact)
# summary(fit4)
# 
# soilmoist <- read.csv("Data//ERG.phytometer.census.csv")
# 
# fit5 <- aovp(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
# summary(fit5)
# t.test(swc ~ Microsite, data=subset(soilmoist, Census=="emergence"))
# 
# micro <- read.csv("Data//ERG.logger.data.csv")
# 
# micro.avg <- micro %>% group_by(Year, Site, Microsite, Rep) %>% summarize(temp=mean(Temp), rh=mean(RH))
# 
# fit6 <- aovp(temp ~ Microsite, data=micro.avg)
# summary(fit6)
# t.test(temp ~ Microsite, data=micro.avg)
# 
# fit7 <- aovp(rh ~ Microsite, data=micro.avg)
# summary(fit7)
# t.test(rh ~ Microsite, data=micro.avg)

Plot level GLMM

m1 <- glmer.nb(Abundance ~ Microsite * gams + (1|Year), data=community.arid)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
car::Anova(m1, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Abundance
##                    Chisq Df Pr(>Chisq)    
## Microsite        16.3613  1  5.234e-05 ***
## gams           1136.3911  1  < 2.2e-16 ***
## Microsite:gams    1.0278  1     0.3107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(m1)
## $Year
##      (Intercept) Micrositeshrub       gams Micrositeshrub:gams
## 2016    13.13187      0.6544052 -0.1298036        -0.005782823
## 2017    11.53434      0.6544052 -0.1298036        -0.005782823
## 
## attr(,"class")
## [1] "coef.mer"
source("rsquared.r")

m2 <- glmer.nb(Richness ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
## Warning in theta.ml(Y, mu, weights = object@resp$weights, limit = limit, :
## iteration limit reached
car::Anova(m2, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Richness
##                            Chisq Df Pr(>Chisq)    
## (Intercept)             140.6109  1  < 2.2e-16 ***
## Microsite                 8.7643  1   0.003072 ** 
## poly(gams, 2)           969.9903  2  < 2.2e-16 ***
## Microsite:poly(gams, 2) 374.2555  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m2)


## richness
p1 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=Richness), color="#181818", size=3, alpha=0.4, width = 0.75,shape = 1, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=Richness), color="#181818", size=2,alpha=0.4, data=subset(community.arid, Microsite=="open"),shape = 2, width = 0.75)  + theme_Publication() +  ylab("species richness") + xlab("Gams index of continentality")+ 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080" , data=subset(community.arid, Microsite=="shrub"), lwd=2) + 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~poly(x,2),aes(x=gams, y=Richness), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+  ylab("Species richness") + xlab("") + annotate("text", x=50,y=6, label="a", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
m3<- lmer(ihs(Biomass) ~ Microsite * poly(gams,2) + (1|Year), data=community.arid)
anova(m3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Microsite                64.53  64.527     1 833.00 210.491 < 2.2e-16 ***
## poly(gams, 2)           511.61 255.805     2 833.43 834.456 < 2.2e-16 ***
## Microsite:poly(gams, 2)  32.08  16.041     2 833.00  52.326 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values
##       R2m       R2c 
## 0.5031296 0.8793048
## biomass
p2 <- ggplot(community.arid) + geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=3, width = 0.75,shape = 1, alpha=0.5, data=subset(community.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=ihs(Biomass)), color="#181818", size=2,shape = 2,alpha=0.5, width = 0.75, data=subset(community.arid, Microsite=="open"))  + theme_Publication() + 
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="shrub"), lwd=2) + 
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=ihs(Biomass)), color="#181818", fill="#80808080", data=subset(community.arid, Microsite=="open"), lwd=2, lty=2)+  ylab("Annual biomass") + xlab("")+ annotate("text", x=50,y=5, label="b", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
m4<- lmer(mpd~ Microsite * gams + (1|Year), data=mpd.gams)
anova(m4)
## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Microsite      0.67216 0.67216     1    24  5.5363 0.02716 *
## gams           0.00313 0.00313     1    24  0.0258 0.87369  
## Microsite:gams 0.54576 0.54576     1    24  4.4952 0.04452 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m4) ## R squared values
##       R2m       R2c 
## 0.2060199 0.2060199
## phylogenetic diversity
p3 <- ggplot(mpd.gams)+ geom_jitter(aes(x=gams, y=mpd), color="#181818", size=3,shape=1, alpha=0.4, data=subset(mpd.gams, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=mpd), color="#181818",alpha=0.4, size=2,shape=2, data=subset(mpd.gams, Microsite=="open"))  +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="shrub"), color="#181818",fill="#80808080", lwd=2)+  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=mpd), data=subset(mpd.gams, Microsite=="open"), color="#181818",fill="#80808080", lwd=2, lty=2)+ ylab("Phylogenetic community dissimilarity") + xlab("Gams index of continentality")+ theme_Publication()+ annotate("text", x=50,y=1.4, label="c", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## native vs non-native
status.gams <- merge(status.data, gams.data, by=c("Year","Site"))

m5<- glmer.nb(native~ Microsite * scale(gams) + (1|Year), data=status.gams) ## need to scale data
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00214781 (tol =
## 0.001, component 1)
car::Anova(m5, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: native
##                         Chisq Df Pr(>Chisq)    
## (Intercept)           565.724  1  < 2.2e-16 ***
## Microsite              23.122  1  1.521e-06 ***
## scale(gams)            67.393  1  2.224e-16 ***
## Microsite:scale(gams)  35.708  1  2.292e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m5)


m6<- glmer.nb(non.native~ Microsite * scale(gams) + (1|Year), data=status.gams)   ## need to scale data
car::Anova(m6, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: non.native
##                          Chisq Df Pr(>Chisq)    
## (Intercept)            12.5323  1   0.000400 ***
## Microsite              13.5445  1   0.000233 ***
## scale(gams)           353.8669  1  < 2.2e-16 ***
## Microsite:scale(gams)   1.4368  1   0.230664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rsquared.glmm(m5)


status.plot <- gather(status.gams, origin, abundance, 5:6)

p4 <- ggplot(subset(status.plot, abundance>0)) + geom_jitter(aes(x=gams, y=abundance), color="#181818", size=3, data=subset(status.plot, Microsite=="shrub"), width = 1, shape=1, alpha=0.3) +geom_jitter(aes(x=gams, y=abundance), width = 1, color="#181818", size=2,alpha=0.3,shape=2, data=subset(status.plot, Microsite=="open"))  + theme_Publication() + 
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance), color="#505050", se=F, data=subset(status.plot, Microsite=="open" & origin=="native"), lwd=2, lty=2)+
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818",  data=subset(status.plot, Microsite=="open" & origin=="non.native"), lwd=2, lty=2)+
 geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#505050", data=subset(status.plot, Microsite=="shrub" & origin=="native"), lwd=2)+
  geom_smooth(method="glm", method.args=list(family="poisson"), formula= y~x,aes(x=gams, y=abundance),se=F, color="#181818",  data=subset(status.plot, Microsite=="shrub" & origin=="non.native"), lwd=2) +  ylab("Annual plant abundance") + xlab("Gams index of continentality")+ annotate("text", x=50,y=310, label="d", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p1,p2,p3,p4)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

### Nutrients 
nutrient.mean <- nutrients.climate %>% group_by(Site, gams, microsite) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean,gams=rep(gams.data[gams.data$Year==2016,"gams"],each=2))

## nitrogen
m5 <- lm(nitrogen ~ poly(gams,2)* microsite, data=nutrient.mean)
summary(m5)
## 
## Call:
## lm(formula = nitrogen ~ poly(gams, 2) * microsite, data = nutrient.mean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2918 -0.9251 -0.2773  0.8012  5.6374 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.595      1.038   1.536 0.163115    
## poly(gams, 2)1                   3.927      3.885   1.011 0.341676    
## poly(gams, 2)2                   3.847      3.885   0.990 0.350969    
## micrositeshrub                   7.676      1.468   5.228 0.000795 ***
## poly(gams, 2)1:micrositeshrub   22.258      5.494   4.052 0.003676 ** 
## poly(gams, 2)2:micrositeshrub   17.833      5.494   3.246 0.011771 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.747 on 8 degrees of freedom
## Multiple R-squared:  0.9298, Adjusted R-squared:  0.8859 
## F-statistic: 21.18 on 5 and 8 DF,  p-value: 0.0002012
anova(m5)
## Analysis of Variance Table
## 
## Response: nitrogen
##                         Df Sum Sq Mean Sq F value    Pr(>F)    
## poly(gams, 2)            2 389.59 194.794  25.817 0.0003239 ***
## microsite                1 206.22 206.223  27.332 0.0007948 ***
## poly(gams, 2):microsite  2 203.35 101.677  13.476 0.0027446 ** 
## Residuals                8  60.36   7.545                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m5) ## R squared values 
##       R2m       R2c 
## 0.8906817 0.8906817
p5 <- ggplot(nutrient.mean) + geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3, shape=16,alpha=0.5, data=subset(nutrient.mean, microsite=="shrub"))+ geom_jitter(aes(x=gams, y=nitrogen), color="#181818", size=3,shape=17, alpha=0.4, data=subset(nutrient.mean, microsite=="open"))+  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080", data=subset(nutrient.mean, microsite=="shrub"), lwd=2) +  
  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=nitrogen), color="#181818", fill="#80808080",data=subset(nutrient.mean, microsite=="open"), lwd=2, lty=2)+ ylab("Nitrogen (ppm)") + xlab("")+ theme_Publication()+ coord_cartesian(ylim=c(-2, 26)) + xlab("")+ annotate("text", x=63,y=25, label="a", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
smc.early <- subset(smc, Census=="emergence")
smc.early <- merge(smc.early, gams.data, by=c("Year","Site"))
smc.early <- smc.early %>% group_by(Year, Site, gams, Microsite) %>% summarize(swc.avg=mean(swc))

## nitrogen
m6 <- lmer(log(swc.avg) ~ poly(gams,2)* Microsite + (1|Year), data=smc.early)
summary(m6)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(swc.avg) ~ poly(gams, 2) * Microsite + (1 | Year)
##    Data: smc.early
## 
## REML criterion at convergence: 16.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75790 -0.67057  0.05989  0.64573  1.52776 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Year     (Intercept) 0.32730  0.5721  
##  Residual             0.09642  0.3105  
## Number of obs: 28, groups:  Year, 2
## 
## Fixed effects:
##                                Estimate Std. Error        df t value
## (Intercept)                    2.183930   0.412962  1.000107   5.288
## poly(gams, 2)1                -4.095881   0.529751 21.483245  -7.732
## poly(gams, 2)2                 0.009355   0.444206 21.039762   0.021
## Micrositeshrub                -0.021937   0.117363 21.000000  -0.187
## poly(gams, 2)1:Micrositeshrub  0.108296   0.621024 21.000000   0.174
## poly(gams, 2)2:Micrositeshrub -0.128193   0.621024 21.000000  -0.206
##                               Pr(>|t|)    
## (Intercept)                      0.119    
## poly(gams, 2)1                1.22e-07 ***
## poly(gams, 2)2                   0.983    
## Micrositeshrub                   0.854    
## poly(gams, 2)1:Micrositeshrub    0.863    
## poly(gams, 2)2:Micrositeshrub    0.838    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(,2)1 pl(,2)2 Mcrsts p(,2)1:
## ply(gms,2)1  0.000                               
## ply(gms,2)2  0.000  0.084                        
## Microstshrb -0.142  0.000   0.000                
## ply(g,2)1:M  0.000 -0.586   0.000   0.000        
## ply(g,2)2:M  0.000  0.000  -0.699   0.000  0.000
anova(m6)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## poly(gams, 2)           8.6913  4.3456     2 21.339 45.0709 2.183e-08 ***
## Microsite               0.0034  0.0034     1 21.000  0.0349    0.8535    
## poly(gams, 2):Microsite 0.0070  0.0035     2 21.000  0.0365    0.9642    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m6) ## R squared values 
##       R2m       R2c 
## 0.5883174 0.9063213
p6 <- ggplot(smc.early) + geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3, shape=16,alpha=0.5, data=subset(smc.early, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=swc.avg), color="#181818", size=3,shape=17,alpha=0.5, data=subset(smc.early, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=swc.avg), data=smc.early, lwd=2, color="black") + ylab("Soil moisture at emergence (%)")  + xlab("")+ annotate("text", x=50,y=35, label="b", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## temperature variability
m7 <- lmer(temp.var~Microsite * poly(gams,2) + (1|Year), data=hobo.arid)
shapiro.test(resid(m7))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m7)
## W = 0.97826, p-value = 0.8069
anova(m7) 
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Microsite                 1.51    1.51     1 21.000  0.2547    0.6190    
## poly(gams, 2)           707.64  353.82     2 21.104 59.7042 2.051e-09 ***
## Microsite:poly(gams, 2)   0.11    0.05     2 21.000  0.0092    0.9908    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m7) ## R squared values 
##       R2m       R2c 
## 0.3627874 0.9560767
p7 <- ggplot(hobo.arid) + geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=16, alpha=0.5, data=subset(hobo.arid, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=temp.var), color="#181818", size=3, shape=17, alpha=0.5, data=subset(hobo.arid, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="open"), lwd=2, lty=2) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=temp.var), color="#181818", fill="#80808080", data=subset(hobo.arid, Microsite=="shrub"), lwd=2)  + ylab("Temperature variation") + xlab("Gams index of continentality")+ annotate("text", x=50,y=125, label="c", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## soil compaction
m8 <- lm(log(compact)  ~  Microsite * gams, data=shrub.clim)
summary(m8)
## 
## Call:
## lm(formula = log(compact) ~ Microsite * gams, data = shrub.clim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37348 -0.12596  0.04006  0.12843  0.33669 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         -0.138777   0.804662  -0.172   0.8665  
## Micrositeshrub      -3.101869   1.137963  -2.726   0.0213 *
## gams                 0.007791   0.010038   0.776   0.4556  
## Micrositeshrub:gams  0.036028   0.014196   2.538   0.0295 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared:  0.6967, Adjusted R-squared:  0.6058 
## F-statistic: 7.658 on 3 and 10 DF,  p-value: 0.006
anova(m8)
## Analysis of Variance Table
## 
## Response: log(compact)
##                Df  Sum Sq Mean Sq F value   Pr(>F)   
## Microsite       1 0.18836 0.18836  3.3176 0.098548 . 
## gams            1 0.75039 0.75039 13.2168 0.004571 **
## Microsite:gams  1 0.36569 0.36569  6.4409 0.029469 * 
## Residuals      10 0.56776 0.05678                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m8) ## R squared values 
##       R2m       R2c 
## 0.6386405 0.6386405
p8 <- ggplot(shrub.clim) + geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=16, alpha=0.5, data=subset(shrub.clim, Microsite=="shrub")) +geom_jitter(aes(x=gams, y=log(compact)), color="#181818", size=3, shape=17, alpha=0.5, data=subset(shrub.clim, Microsite=="open"))  + theme_Publication()+ stat_smooth(method="lm", formula= y~x, aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="open"), lwd=2, lty=2) +  stat_smooth(method="lm", formula= y~x,aes(x=gams, y=log(compact)), color="#181818", fill="#80808080", data=subset(shrub.clim, Microsite=="shrub"), lwd=2)  + ylab("Soil compaction") + xlab("Gams index of continentality")+ annotate("text", x=63,y=1.1, label="d", size=8) 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p5,p6,p7,p8)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

Phytometer with GAMS

census <- read.csv("Data/ERG.phytometer.census.csv")
census[is.na(census)] <- 0
census.arid <- merge(census,gams.data, by=c("Year","Site"))

## treatment plots
census.long <- gather(census.arid, species, biomass, 12:14)
census.long[,"species"] <- ifelse(census.long[,"species"] =="Phacelia.biomass", "P. tanacetifolia",
                                  ifelse(census.long[,"species"]=="Plantago.biomass","P. insularis","S. columbariae"))
census.long[,"species"] <- as.factor(census.long$species)

## microsite average
census.micro <- census.long %>% filter(biomass>0) %>%  group_by(Microsite, species) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

p9 <- ggplot(census.micro, aes(x=species, y=Biomass, fill=Microsite)) +   geom_bar(position=position_dodge(), stat="identity", colour="black") + scale_fill_brewer(palette="Greys")+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication() + theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## nutrient average
census.nut <- census.long %>% filter(biomass>0) %>%  group_by(Nutrient, species) %>%  summarize(Biomass=mean(biomass), error=se(biomass))

p10 <- ggplot(census.nut, aes(x=species, y=Biomass, fill=Nutrient)) +   geom_bar(position=position_dodge(), stat="identity",colour="black") + scale_fill_brewer(palette="Greys")+  geom_errorbar(aes(ymin=Biomass-error, ymax=Biomass+error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ theme_Publication()+ theme(axis.text.x = element_text(size=16, face="italic"),axis.text.y = element_text(size=16),axis.title = element_text(size=20),legend.text=element_text(size=18))
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
## logistic first
census.long[,"presence"] <- ifelse(census.long$biomass>0,1,0)
## Drop panoche 2017 because of cold
census.npan <- census.long %>% filter(Site!="PanocheHills" | Year!="2017")


m1.l <- glmer(presence~Microsite * gams * Nutrient + (1|Year),data=subset(census.npan, species=="P. tanacetifolia"), family=binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00301814 (tol =
## 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
car::Anova(m1.l, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                           Chisq Df Pr(>Chisq)    
## Microsite               23.7693  1  1.086e-06 ***
## gams                     6.7832  1   0.009202 ** 
## Nutrient                 0.0047  1   0.945204    
## Microsite:gams           3.2156  1   0.072938 .  
## Microsite:Nutrient       1.6347  1   0.201055    
## gams:Nutrient            0.0000  1   0.995620    
## Microsite:gams:Nutrient  0.0062  1   0.937144    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.long, species=="P. tanacetifolia" & presence==1))
shapiro.test(resid(m1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m1)
## W = 0.99225, p-value = 0.2595
anova(m1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq Mean Sq NumDF DenDF F value
## Microsite                         22.121  22.121     1   221 12.2132
## poly(gams, 2)                    157.099  78.549     2   221 43.3685
## Nutrient                           5.710   5.710     1   221  3.1524
## Microsite:poly(gams, 2)           15.688   7.844     2   221  4.3307
## Microsite:Nutrient                 0.078   0.078     1   221  0.0429
## poly(gams, 2):Nutrient             2.854   1.427     2   221  0.7878
## Microsite:poly(gams, 2):Nutrient   1.406   0.703     2   221  0.3883
##                                     Pr(>F)    
## Microsite                        0.0005733 ***
## poly(gams, 2)                    < 2.2e-16 ***
## Nutrient                         0.0771914 .  
## Microsite:poly(gams, 2)          0.0142925 *  
## Microsite:Nutrient               0.8361223    
## poly(gams, 2):Nutrient           0.4561215    
## Microsite:poly(gams, 2):Nutrient 0.6786980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m1) ## R squared values 
##      R2m      R2c 
## 0.346015 0.346015
m2.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="P. insularis"), family=binomial)
car::Anova(m2.l, type=2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: presence
##                                   Chisq Df Pr(>Chisq)  
## Microsite                        1.4664  1    0.22591  
## poly(gams, 2)                    8.3747  2    0.01519 *
## Nutrient                         0.2723  1    0.60178  
## Microsite:poly(gams, 2)          4.8534  2    0.08833 .
## Microsite:Nutrient               1.6652  1    0.19691  
## poly(gams, 2):Nutrient           2.9034  2    0.23418  
## Microsite:poly(gams, 2):Nutrient 0.4909  2    0.78236  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="P. insularis" & presence==1))
shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.99193, p-value = 0.5364
anova(m2)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                  Sum Sq Mean Sq NumDF  DenDF F value
## Microsite                         0.569  0.5693     1 141.23  0.5147
## poly(gams, 2)                    36.153 18.0765     2 137.31 16.3429
## Nutrient                          8.958  8.9576     1 141.27  8.0986
## Microsite:poly(gams, 2)           0.858  0.4289     2 141.51  0.3878
## Microsite:Nutrient                0.110  0.1096     1 141.02  0.0991
## poly(gams, 2):Nutrient            2.373  1.1866     2 141.20  1.0728
## Microsite:poly(gams, 2):Nutrient  0.654  0.3272     2 141.03  0.2958
##                                     Pr(>F)    
## Microsite                          0.47430    
## poly(gams, 2)                    4.297e-07 ***
## Nutrient                           0.00509 ** 
## Microsite:poly(gams, 2)            0.67929    
## Microsite:Nutrient                 0.75338    
## poly(gams, 2):Nutrient             0.34481    
## Microsite:poly(gams, 2):Nutrient   0.74439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m2) ## R squared values 
##       R2m       R2c 
## 0.2802814 0.5695941
m3.l <- glmer(presence~Microsite * poly(gams,2) * Nutrient + (1|Year),data=subset(census.npan, species=="S. columbariae"), family=binomial)
car::Anova(m3.l, type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: presence
##                                    Chisq Df Pr(>Chisq)    
## (Intercept)                      55.9936  1  7.271e-14 ***
## Microsite                         1.2916  1     0.2558    
## poly(gams, 2)                    18.5593  2  9.330e-05 ***
## Nutrient                          0.0728  1     0.7873    
## Microsite:poly(gams, 2)           1.7247  2     0.4222    
## Microsite:Nutrient                2.3865  1     0.1224    
## poly(gams, 2):Nutrient            0.8524  2     0.6530    
## Microsite:poly(gams, 2):Nutrient  0.5872  2     0.7456    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- lmer(log(biomass)~Microsite * poly(gams,2) * Nutrient + (1|Year), data=subset(census.long, species=="S. columbariae" & biomass>0))
shapiro.test(resid(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m3)
## W = 0.99363, p-value = 0.1532
anova(m3)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                  Sum Sq Mean Sq NumDF  DenDF F value
## Microsite                         0.350  0.3503     1 333.03  0.2070
## poly(gams, 2)                    61.028 30.5139     2 327.98 18.0272
## Nutrient                         18.797 18.7973     1 333.12 11.1052
## Microsite:poly(gams, 2)           1.304  0.6521     2 333.03  0.3853
## Microsite:Nutrient                0.386  0.3859     1 333.05  0.2280
## poly(gams, 2):Nutrient            2.914  1.4569     2 333.05  0.8607
## Microsite:poly(gams, 2):Nutrient  3.356  1.6780     2 333.01  0.9913
##                                     Pr(>F)    
## Microsite                        0.6494522    
## poly(gams, 2)                    3.733e-08 ***
## Nutrient                         0.0009576 ***
## Microsite:poly(gams, 2)          0.6805772    
## Microsite:Nutrient               0.6333375    
## poly(gams, 2):Nutrient           0.4237872    
## Microsite:poly(gams, 2):Nutrient 0.3721737    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r.squaredGLMM(m3) ## R squared values 
##       R2m       R2c 
## 0.1506253 0.2337410
## summary presence
mean.occ <- census.npan %>% group_by(gams,species, Microsite) %>% summarize(presence=mean(presence))


p11 <- ggplot(census.npan, aes(x=gams, y=presence, fill=species, color=species))  + theme_Publication()+  
  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2",  data=subset(census.npan, species=="P. tanacetifolia" & Microsite=="shrub"), lwd=2)+  
  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~x,aes(x=gams, y=presence),fill="#0daec2", color="#0daec2",  data=subset(census.npan, species=="P. tanacetifolia"& Microsite=="open"), lwd=2, lty=2) +  geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence), color="#fb6f70", fill="#fb6f7050",  data=subset(census.npan, species=="P. insularis"), lwd=2)+ geom_smooth(method="glm", method.args=list(family="binomial"), formula= y~poly(x,2),aes(x=gams, y=presence),fill="#b7d73850", color="#b7d738",  data=subset(census.npan, species=="S. columbariae"), lwd=2)+ ylab("Probability of occurrence") + xlab("Gams index of continentality")
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
p12 <- ggplot(census.long, aes(x=gams, y=biomass, fill=species, color=species))  + theme_Publication() +  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)), color="#fb6f70", fill="#fb6f7050",  data=subset(census.long, species=="P. insularis" & biomass>0), lwd=2, fullrange=TRUE) + geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#b7d73850", color="#b7d738",  data=subset(census.long, species=="S. columbariae"& biomass>0), lwd=3)+
  geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2",  data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="shrub"), lwd=2, fullrange=TRUE)+geom_smooth(method="lm", formula= y~poly(x,2),aes(x=gams, y=log(biomass)),fill="#0daec2", color="#0daec2",  data=subset(census.long, species=="P. tanacetifolia" & biomass>0 & Microsite=="open"), lwd=2,lty=2, fullrange=TRUE)+ ylab("Phytometer biomass") + xlab("Gams index of continentality") 
## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing
grid.arrange(p9,p11,p10, p12)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database